Programmer Power Tools/Programmer Power Tools.iso/microcrn/issue_40.arc/DAIMS.ARC/MATRIX2.CPP

C/C++ Source or Header | 1988-02-10 | 13.0 KB | 501 lines

#include "matrix.hxx"
#include <math.h>
#include <ctype.h>
#include <time.h>
#ifndef DOS
#include <sys/param.h>
#endif


/*
-*++ matrix::set_row(): sets a row to a value
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

void matrix::set_row(int row, double value)
{
  for(int i = 0; i < cols(); i++)
    (*this)[row][i] = value;
}


/*
-*++ matrix::set_column(): sets a column to a value
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

void matrix::set_column(int column, double value)
{
  for(int i = 0; i < cols(); i++)
    (*this)[i][column] = value;
}


/*
-*++ matrix::switch_rows(): switches two rows of a matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** Could be optimized.
** ++*)
**
** (*++ detailed:
** ++*)
*/

void matrix::switch_rows(int row1, int row2)
{
  matrix temp(1,cols());
  for(int col = 0; col < cols(); col++)
    temp[0][col] = (*this)[row1][col]; /* temporarily store row 1 */
  for(col = 0; col < cols(); col++)
    (*this)[row1][col] = (*this)[row2][col]; /* move row2 to row1 */
  for(col = 0; col < cols(); col++)
    (*this)[row2][col] = temp[0][col]; /* move temp to row2 */
}


/*
-*++ matrix::switch_columns(): switches two columns of a matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

void matrix::switch_columns(int col1, int col2)
{
  matrix temp(rows(),1);
  for(int row = 0; row < rows(); row++)
    temp[row][0] = (*this)[row][col1]; /* temporarily store col 1 */
  for(row = 0; row < rows(); row++)
    (*this)[row][col1] = (*this)[row][col2]; /* move col2 to col1 */
  for(row = 0; row < rows(); row++)
    (*this)[row][col2] = temp[row][0]; /* move temp to col2 */
}


/*
-*++ matrix::row_multiply_add(): adds a multiple of one row to another
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

void matrix::row_multiply_add(double multiplier, int reference_row, int changing_row)
{
  for(int col = 0; col < cols(); col++)
    (*this)[changing_row][col] += (*this)[reference_row][col] * multiplier;
}


/*
-*++ matrix::col_multiply_add(): adds a multiple of one col to another
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

void matrix::col_multiply_add(double multiplier, int reference_col, int changing_col)
{
  for(int row = 0; row < rows(); row++)
    (*this)[row][changing_col] += (*this)[row][reference_col] * multiplier;
}


/*
-*++ matrix::min(): returns the minimum element in the matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

double matrix::min()
{
  double temp;
  if(rows() <= 0 || cols() <= 0) error("bad matrix size for min()");
  double minimum = (*this)[0][0];
  for (int row = 0; row < rows(); row++)
    for(int col = 0; col < cols(); col++)
      if ((temp = (*this)[row][col]) < minimum)
    minimum = temp;
  return minimum;
}


/*
-*++ matrix::max(): returns the maximum element in the matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

double matrix::max()
{
  double temp;
  if(rows() <= 0 || cols() <= 0) error("bad matrix size for max()");
  double maximum = (*this)[0][0];
  for (int row = 0; row < rows(); row++)
    for(int col = 0; col < cols(); col++){
      if ((temp = (*this)[row][col]) > maximum)
    maximum = temp;
    }
  return maximum;
}


/*
-*++ matrix::mean(): returns the mean of all the matrix elements
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

double matrix::mean()
{
  double sum = 0;
  for (int row = 0; row < rows(); row++)
    for(int col = 0; col < cols(); col++)
      sum += (*this)[row][col];
  return sum/(row * col);
}


/*
-*++ matrix::variance(): returns the statistical variance of all elements
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

double matrix::variance()
{
  double s_squared = 0;
  double mn = mean();
  for (int row = 0; row < rows(); row++)
    for(int col = 0; col < cols(); col++){
      double temp = (*this)[row][col] - mn;
      temp *= temp;
      s_squared += temp;
    }
  s_squared /= row * col -1;    /* number of elements minus one */
  return s_squared;
}


/*
-*++ matrix::transpose(): returns the transpose matrix (must be square)
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix &  matrix::transpose()
{
  if(rows() != cols()) error("matrix must be square to transpose!\n");
  matrix & trans= *new matrix(rows(),cols());
  for (int row = 0; row < rows(); row++)
    for(int col = 0; col < cols(); col++)
      trans[col][row] = (*this)[row][col];
  return trans;
}


/*
-*++ matrix::scale(): scale a matrix (used in LU decomposition)
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix & matrix::scale()
{
  double temp;
  if(rows() <= 0 || cols() <= 0) error("bad matrix size for scale()");
  if(rows() != cols()) error("matrix must be square for scale()");
  matrix & scale_vector= *new matrix(rows());
  for (int col = 0; col < cols(); col++){
    double maximum = 0;
    for(int row = 0; row < rows(); row++)
      if ((temp = fabs((*this)[row][col])) > maximum)
    maximum = temp;        /* find max column magnitude in this row */
    if(maximum == 0) error("singular matrix in scale()");
    scale_vector[col][0] = 1/maximum; /* save the scaling */
  }
  return scale_vector;
}


/*
-*++ matrix::bitcopy(): make an image of a matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

void matrix::bitcopy(matrix& from, matrix& to)
{
  if(from.rows() != to.rows() || from.cols() != to.cols())
    error("matrices must be the same dimensions for bitcopy()");
  for(int row = 0; row < from.rows(); row++)
    for(int col = 0; col < from.cols(); col++)
      to[row][col] = from[row][col];
}


/*
-*++ matrix::lu_decompose(): returns LU decomposition matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: indx is an output vector which records the row permutation
** effected by the partial pivoting, d is output as +-1 depending on whether
** the number of row interchanges was even or odd, respectively.  This routine
** is used in combination with lu_back_subst to solve linear equations or
** invert a matrix.
**     Not much in the way of comments here because I don't really know what
** I'm doing.  I'm following "Numerical Recipes: The Art of Scientific
** Computing," by Press and Flannery, Cambridge Press 1986 pp. 31-39.
** ++*)
*/

matrix & matrix::lu_decompose(matrix& indx, int& d )
{
  if(rows() != cols()) error("Matrix must be square to L-U decompose!\n");
  d = 1;            /* parity check */
  int row,col,k,col_max; /* counters */
  double dum;            /* from the book -- I don't know significance */
  double sum;
  double maximum;
  matrix & lu_decomp = *new matrix(rows(), cols());
  bitcopy(*this,lu_decomp);    /* direct copy of the original matrix */
//  lu_decomp = *this + 0;    /* makes a copy of this (gets around ref)  */
  matrix scale_vector(lu_decomp.scale()); /* scale the matrix */

  for(row = 0; row < rows(); row++){ /* The loop over columns of Crout's method */
    if (row > 0) {
      for (col = 0; col < row; col++) { /* eqn 2.3.12 except for row=col */
    sum = lu_decomp[row][col];
    if(col > 0) {
      for(k = 0; k < col; k++)
        sum -= lu_decomp[row][k] * lu_decomp[k][col];
      lu_decomp[row][col] = sum;
    }
      }
    }
    maximum = 0; /* Initialize for the search for the largest pivot element */
    for(col=row; col < cols(); col++){ /* i=j of eq 2.3.12 & i=j+1..N of 2.3.13 */
      sum = lu_decomp[row][col];
      if(row > 0){
    for(k=0; k < row; k++)
      sum -=  lu_decomp[k][col] * lu_decomp[row][k];
    lu_decomp[row][col] = sum;
      }
      dum = scale_vector[col][0] * fabs(sum);/*figure of merit for pivot*/
      if (dum >= maximum){ /* is it better than the best so far? */
    col_max = col;
    maximum = dum;
      }
    }
    col--; /* fix -- maybe difference between fortran and C for-loops. */
       /* I think maybe Fortran leaves the counter at the high */
       /* value while C increments it one past the high value */
    if(row != col_max) {       /* Do we need to interchange rows? */
      lu_decomp.switch_columns(col_max,row); /* Yes, do so... */
      d *= -1;               /* ... and change the parity of d */
      dum = scale_vector[col_max][0]; /* also interchange */
      scale_vector[col_max][0] = scale_vector[col][0]; /* the scale factor */
      scale_vector[row][0] = dum; /* Keffer did it this way, not P&F */
    }
    indx[row][0] = col_max;
    if(row != rows() -1){  /* Now, finally, divide by the pivot element */
      if(lu_decomp[row][row] == 0) lu_decomp[row][row] = TINY;
      /* If the pivot element is zero the matrix is singular (at least to the
     precision of the algorithm).  For some applications on singular matrices,
     it is desirable to substitute TINY for zero/ */
      dum = 1/lu_decomp[row][row];
      for(col=row+1; col < cols(); col++)
    lu_decomp[row][col] *= dum;
    }
  }
  if(lu_decomp[rows()-1][ cols()-1] == 0)
    lu_decomp[rows()-1][cols()-1] = TINY;
  return lu_decomp;
}



/*
-*++ matrix::lu_back_subst(): Solves the set of N linear equations A*X = B
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: Here "this" is the LU-decomposition of the matrix A,
** determined by the routine lu_decompose(). Indx is input as the permutation
** vector returned by lu_decompose().  B is input as the right-hand side
** vector B, and returns with the solution vector X.  This routine takes into
** account the possibility that B will begin with many zero elements, so it is
** efficient for use in matrix inversion.   See pp 36-37, Press and Flannery.
** ++*)
*/

    void matrix::lu_back_subst(matrix& indx, matrix& b)
{
  if(rows() != cols())
    error ("non-square lu decomposition matrix in lu_back_subst()");
  if(rows() != b.rows())
    error("wrong size B vector passed to lu_back_subst()");
  if(rows() != indx.rows())
    error("wrong size indx vector passed to lu_back_subst()");
  int row,col,ll;
  int ii = 0;
  double sum;
  for(col=0;col < cols(); col++){
    ll= (int)indx[col][0];
    sum = b[ll][0];
    b[ll][0] = b[col][0];
    if (ii >= 0)
      for(row = ii; row < col; row++)
    sum -= (*this)[row][col] * b[row][0];
    else if(sum != 0)
      ii = col;
    b[col][0] = sum;
  }
  for(col = cols() -1; col >= 0; col--){
    sum = b[col][0];
    if (col < cols() -1)
      for (row = col + 1; row < rows(); row++)
    sum -= (*this)[row][col] * b[row][0];
    b[col][0] = sum/(*this)[col][col]; /* store a component of the soln vector X*/
  }
}


/*
-*++ matrix::copy_column(): copy one column onto another
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed: copy the from_col of m to the to_col of this
** ++*)
*/

void matrix::copy_column(matrix& m, int from_col, int to_col)
{
  if(rows() != m.rows()) error("number of rows must be equal for copy_column()");
  for(int row=0; row < rows(); row++)
    (*this)[row][to_col] = m[row][from_col];
}



/*
-*++ matrix::inverse(): returns the inverse matrix of a square matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

matrix & matrix::inverse()
{
    if(rows() != cols()) error("matrix must be square for inverse()");
    matrix Y("I",rows());        /* create an identity matrix */
    matrix indx(cols());        /* create the "index vector" */
    matrix B(cols());        /* see pp 38. */
    int d;
    matrix & decomp = *new matrix(lu_decompose(indx,d));
                /* perform the decomposition once */
    for(int col = 0; col < cols(); col++){
    B.copy_column(Y,col,0);
    decomp.lu_back_subst(indx,B);
    Y.copy_column(B,0,col);
    }
    return Y.transpose();
}


/*
-*++ matrix::determinant(): returns the determinant of a square matrix
**
** (*++ history:
**      4 Dec 87    Bruce Eckel    Creation date
** ++*)
**
** (*++ detailed:
** ++*)
*/

double matrix::determinant()
{
  if(rows() != cols()) error("matrix must be square for determinant()");
  matrix indx(cols());        /* create the "index vector" */
  matrix B(cols());        /* see pp 38. */
  int d;
  matrix decomp(lu_decompose(indx,d)); /* perform the decomposition once */
  double determinant = d;
  for(int i=0; i < cols() ; i++)
    determinant *= decomp[i][i];
  return determinant;
}